Explosive Percolation is Continuous, but with Unusual Finite Size Behavior 
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We study four Achlioptas type processes with "explosive" percolation transitions. All transi- 
tions are clearly continuous, but their finite size scaling functions are not entire holomorphic. The 
distributions of the order parameter, the relative size Smax/iV of the largest cluster, are double- 
humped. But - in contrast to first order phase transitions - the distance between the two peaks 
decreases with system size N as N~ v with n > 0. We find different positive values of /3 (defined 
via (s max /7V) ~ (p — Pc) 13 for infinite systems) for each model, showing that they are all in different 
universality classes. In contrast, the exponent (defined such that observables are homogeneous 
functions of (p — p c )N e ) is close to - or even equal to - 1/2 for all models. 

PACS numbers: 64.60.ah, 05.70.Jk, 89.75.Da, 05.40.-a 



Percolation is a pervasive concept in statistical physics 
and probability theory and has been studied in extenso 
in the past. It came thus as a surprise to many, when 
Achlioptas et al. [l| claimed that a seemingly mild modi- 
fication of standard percolation models leads to a discon- 
tinuous phase transition - named "explosive percolation" 
(EP) by them - in contrast to the continuous phase tran- 
sition seen in ordinary percolation. Following [lj there 
appeared a flood of papers (3-[20j studying various as- 
pects and generalizations of EP. In all cases, with one ex- 
ception [20( , the authors agreed that the transition is dis- 
continuous: the "order parameter" , defined as the frac- 
tion of vertices/sites in the largest cluster, makes a dis- 
crete jump at the percolation transition. In the present 
paper we join the dissenting minority and add further 
convincing evidence that the EP transition is continuous 
in all models, but with unusual finite size behavior. 

From the physical point of view, the model seems some- 
what unnatural, since it involves non-local control (there 
is a 'supervisor' who has to compare distant pairs of 
nodes to chose the actual bonds to be established JUJ). 
Also, notwithstanding no realistic applications have 
been proposed. It is well known that the usual concept of 
universality classes in critical phenomena is invalidated 
by the presence of long range interactions. Thus it is not 
surprising that a percolation model wi th g lobal control 
can show completely different behavior [22( . 

Usually, e.g. in thermal equilibrium systems, discon- 
tinuous phase transitions are identified with "first order" 
transitions, while continuous transitions are called "sec- 
ond order" . This notation is also often applied to per- 
colative transitions. But EP lacks most attributes - ex- 
cept possibly for the discontinuous order parameter jump 
- considered essential for first order transitions. None 
of these other attributes (cooperativity, phase coexis- 
tence, and nucleation) is observed in Achlioptas type pro- 
cesses, although they are observed in other percolation- 
type transitions [jjj ]. Thus EP should never have been 
viewed as a first order transition, and it is gratifying that 



it is also not discontinuous. 

Apart from the behavior of the average value (m) of the 
order parameter m, phase transitions can also be char- 
acterized by the distribution P Pi at(to) of m in finite sys- 
tems, where p is the control parameter and N measures 
the system size. For infinite N, (m) jumps at p = p c 
if the transition is discontinuous, while it varies contin- 
uously with a power law singularity (m) <~ (p — Pc) 13 for 
a continuous transition. The distribution P t 
criticality scales, for continuous transitions, as 




P, 



<N {m) ~ TfiftmN*), 



(1) 



where n = j3/{dv) for standard thermal second order 
phase transitions. The universal function f(z) might be 
double- humped, as in the Ising model [24| . But then, 
as N — > oo, the dip between the humps usually does not 
deepen and the horizontal distance between them shrinks 
to zero so that Pp= Pc .N(m) becomes single- humped. 

Equation (1) is directly related to the finite size scaling 
(FSS) of (m) Q, 



(to) ~ (p - p c fg[{p - Pc)N 



(2) 



where the universal scaling function g(z) is analytic at 
all finite z, reflecting the fact that the critical point was 
the only singularity of the partition function, before it 
was re gula rized by Eq.©. Notice that the usual FSS 
ansatz [25| involves the linear system size L instead of N 
with = l/(du), where d is the dimension and v is the 
correlation length exponent. 

In typical first order transitions, in contrast, 
Pp= Pa ,N(, m ) is double- humped with a deepening valley 
between the two peaks. The distance between the peaks 
tends to a positive constant which is equal to the jump in 
(m) . The depth of the valley between the peaks reflects 
the fact that values of m between the peaks correspond 
to systems with two co-existing phases and an interface 
between them that costs energy and is disfavored. As 
a consequence, systems with first order transitions typi- 
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FIG. 1. (Color online) Distributions of the order parameter 
Smax/N for four EP models. They are shown at the effective 
critical point, defined such that both peaks have the same 
height. Normalization is such that their height is I. For the 
largest systems, curves were approximated by cubic splines to 
make them smooth. 



cally do not show FSS (unless the interface energy does 
not increase with system size pfjj ) . 

In percolation, usually the relative size of the largest 
cluster, m = s max /N, is taken as an order parameter. 
Here, N is the number of nodes, and s max /N — > for 
p < p c and N — > oo. In 0, [l6| it was observed that 
P p=p „ Ar(m) is strongly double-peaked in EP transitions. 
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In 16[ this was also backed by careful measurements of 



TABLE I. Critical points and critical exponents for the four 
models. The ©i are different estimates of the exponent 0: <di 
is obtained from the scaling relation = i]+//3, O2 is obtained 
from a data collapse in the slightly supercritical region where 
(m) ~ m+, and C onj is the conjectured exact value. For the 
da Costa model, p c is taken from [2(J. For the other models 
it is obtained from plots analogous to the inset in Fig. 4. 



the depth of the valley between the peaks, which indeed 
lowered with increasing N. This was taken as a clear in- 
dication for the transition being first order and for phase 
coexistence. Notice that the latter is not justified since 
Smax/JV is, in contrast to the local order parameters in 
thermal systems, a global quantity and cannot be used 
to characterize any part of a large system. Rather, the 
structure of P p = Pa ,N{m) in EP reflects the suddenness 
of the transition, combined with a scatter of the precise 
p- values where individual systems acquire giant clusters. 
At p- values where both peaks have the same height, it is 
much more likely to find either no giant cluster or a fully 
developed one, than to find a half-grown giant cluster. 
Hence, the two peaks are more reminiscent of systems 
without self-averaging [22| than of phase coexistence. 

While the two peaks prove the suddenness of the tran- 
sition that was claimed as a hallmark of EP, they do 
not yet prove that EP is discontinuous. For that, one 
must also show that the distance between the peaks does 
not vanish for N — > 00. In order to check this, we have 
made extensive simulations of four models: The original 
product rule of [H], denoted in the following as "PR"; 
The product rule on 2-d square lattices 0, H| with helical 
boundary conditions ( "2d" ) ; The 'adjacent edge' rule Q 
("AE"); And the rule of |o| ("da Costa"). For more de- 
tails on the simulations, see the supplementary material 
(SM). 

Distributions P p ,N(jn) for these models are shown in 
Fig. 1. In all cases p was chosen such that both peaks 
have equal height (set arbitrarily to 1). The extrapo- 
lations of these values for N — » 00 are given in Table 
1. They agree within errors with the critical p c values 
quoted in the literature. We see that in each case the 
valley between the peaks deepens with increasing N [l(| , 
but at the same time both peaks shift to the left. Among 
the three off-lattice models, the AE model (the least non- 
local) shows the fastest peak shifting and slowest valley 
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FIG. 2. (Color online) Data collapses for the two peaks in the 
order parameter distribution for the da Costa model. Colors 
and line styles are the same as in Fig. 1. 



deepening, while the opposite is true for the da Costa 
model. In all cases this shift is compatible with power 
laws 



m± 



N -v± 



(3) 



where m+ (to_) is the position of the right (left) peak 
at the critical point. In all cases < ?/+ < r/_ (sec Table 
1), i.e. the right peak moves slower than the left one. 
Therefore the distance between the peaks increases for 
small N, but has finally to decrease ~ N~ 7, + . Since this 
distance is asymptotically proportional to the maximum 
of the variance of to [28| , we find that the variances first 
increase with N (in agreement with Q), but ultimately 
must decrease. 

As shown in Fig. 2 for the da Costa model, not only the 
positions of the peaks scale, but also their widths. This 
indicates that the asymptotic scenario is two well sepa- 
rated peaks with TV— independent shapes whose widths 
are proportional to their positions. If wc switch from 
defining p c by equal peak heights to equal peak areas 
[28] and allow weak convergence for TV — > oo (in con- 
trast to the usual assumption of pointwise convergence; 
see SM) the full distributions at p c (iV) then show asymp- 
totic scaling 

p Pe mA™)~ NV+ f( mNV+ ) ( 4 ) 

with the scaling function f(x) consisting of a finite width 
right hand peak and a 8— peak at x = 0. 

For p strictly larger than p c (N), only the right hand 
peak dominates the average (to). We then expect only 
small finite size scaling corrections to its asymptotic val- 
ues, i.e. we expect the curves (to) p jv for different N 
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FIG. 3. (Color online) Log-log plot of the average order pa- 
rameter for the da Costa model versus p — p c , for six different 
values of N. One sees clearly a common part with slope ft 
(indicated also by the straight line), from which curves for 
different N deviate later and later, as N increases. The in- 
set shows the collapse of these data as predicted by Eq. |[2j. 
While <d is fitted, both p and p c are taken from 



to coincide for p > p c {N) on a common curve (m) p . 
Since the scenario in this regime is not much different 
from other critical phenomena this should be a power law 
(to) ~ (p — Pc)@ that holds in the range m + < (to) <C 1. 
Measured values of /3 are given in Table 1. For the da 
Costa model the agreement with is perfect. Assum- 
ing Eq. (121), it follows that 8 = ?y+//3. Values of 6 ob- 
tained from this, denoted as 0i, are slightly smaller than 
1/2 for all models (see Table 1). 

Deviations from this common power law are expected 
to set in when (to) decreases below m + . The data for 
the da Costa model are shown in Figs. 3. For all p > 
p c (except for very small values of z = (p — p c )N ), 
these deviations arc fully described by the FSS ansatz in 
Eq. In Figs. 3 we chose O so that the collapse is best 
at (to) ss to_|_, resulting in the value O2 quoted in Table 1. 
For the other models the data collapse is similarly good, 
except for the 2d model where it is worse (see SM). For 
all models, 02 is slightly larger than 0i. 

The fact that f(z) in Eq. ((4]) contains a <5-peak at its 
leftmost extremity z = implies that g(z) in Eq.([2|) must 
vanish for all z below some value zq < 0, which in turn 
means that g(z) must have a singularity at zq. Indeed, 
Fig. 4 shows that the values of g(z) for z < — 1 approach 
very fast with increasing N, implying — 1 < z < 
(the latter is also true for the other models). We cannot 
exclude the possibility the curves in Fig. 4 approach a 
pure power law az@ (dashed red line) in the limit N — > 00. 

The blow-up of the region around z = shown in the 
inset in Fig. 4 hints at a power law (to)| p=Pc ~ N~ m with 
770 = 0.0356(8) > rj + (see also SM). The same is qualita- 
tively true for the other models, where always 770 > rj+ 
(see Table 1). We see therefore that z = is no longer in 
the realm of uniform pointwise convergence to the FSS 
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FIG. 4. (Color online) Doubly linear plot of the same data 
shown in Fig. 3, but extended to values p < p c . Here we used 
<d — 1/2, which gives worse data collapse for p > p c , but 
vastly more systematic behavior for p < p c . The inset shows 
a blow-up of the region around p = p c , with logarithmic y- 
axis. The decrease of the curves at z = with N suggests 
that zq = 0, and that a new power law holds for p — p c . 



ansatz, and that therefore zq = 0. We should finally men- 
tion that we used 8 = 0.5 in Fig. 4, a value in between 
0! and 8 2 , as it gives the most systematic behavior for 
z < 0. The same is true the other off-lattice models (but 
not for the 2d model, see SM), whence we conjecture that 
8 = 0.5 for them. 

The singularity of g(z) at z = implies also that 
one cannot expect the effective critical points to scale 
as p c (N) — p ~ N~®. Results obtained for the da Costa 
model, with p c (N) defined via equal peak masses, are 
shown in the SM. They indicate that p c (N) - p ~ N~ s 
with 5 = 0.9(1) > 8. The agreement with the predic- 
tion 5 = 0.818(1) of - based on "standard scaling 
relations" - seems fortitious. 

In this paper we do not present a detailed theory for 
the convergence to g(z) for z < 0, in particular we do not 
explain how t/q and 5 are related to the exponent rj- . It 
could be that such a theory can be formulated more easily 
using cither (logs max ) or (l/s max ) as an order parameter. 
But this would be beyond the scope of the present paper. 

In summary, we have shown that at least four models of 
explosive percolation, including the original product rule 
of Achlioptas et at [l|, have continuous transitions. Each 
is in a different universality class, but all of them show 
unusual finite size behavior with a non-analytic scaling 
function. They all show double-peaked order parameter 
distributions with the sharpness of the peaks increasing 
with system size, and different scaling laws for the width 
of the scaling region ( <~ N~ & ) and for the shift of the 
effective p c (N). It would be interesting to see whether 
similar scaling holds in other percolation models with 
supposedly discontinuous transitions that are not explic- 
itly related to Achlioptas- like dynamics 0, Q, H^. It 
could be that the features found in the present paper 



arise from the specific non-locality of the Achlioptas pro- 
cess, and that this is why it was not seen previously in 
other critical phenomena. 

We are indebted to Bob Ziff and Liang Tian for most 
useful correspondence. 
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SIMULATION DETAILS 

All simulations reported in the paper were made using modified versions of the fast Newman-Ziff algorithm 
on a Linux workstation cluster. System sizes varied between N = 2 10 and N = 2 26 (w 6.7 x 10 7 ) (N is the number 
of nodes). For the smaller systems « 10 8 realizations were made for each model, and for the largest systems this 
number was still > 10 4 . The control parameter p is defined as in the references where the 4 models were introduced, 
as p = L/N where L is the number of links. 

Data were actually collected for fixed n, where n is the number of clusters - more precisely, in order to reduce the 
data files, n was binned (typically with An = 1 for smallest N and An = 2 8 for largest N). The values of p quoted in 
the paper are average values over these bins. Since most clusters in all four models are trees, except when p is very 
large, there are very small fluctuations of p for fixed n, and (p) depends smoothly on n. Moreover, test runs showed 
that the dependence of s max on n is at least as crisp as the dependence on p, i.e. n is actually the more relevant 
control parameter. 

Mass distributions (Figs. 1 and 2 in the main paper) are obtained by binning, with typically 200 to 500 bins, and 
with bin sizes slowly increasing with s max in order to take into account that the left hand peaks in Fig. 1 are sharper 
than the right hand peaks. For small N the distributions shown in Figs. 1 and 2 are the raw data, modified just by 
interpolating between neighboring n bins to obtain exactly equally high peaks (usually, no bin will have two peaks 
which have exactly equal height; by "interpolating" we mean taking weighted linear averages of the two histograms) 
and by normalizing them. For the largest N this would have given too noisy plots, and cubic splines were used to 
make the plots more smooth. 

Unless otherwise noted, p c values are those in Table 1. They were determined by having best power laws ~ N~ m 
for (m) at p = p c . 



MODIFIED FINITE SIZE SCALING 

The finite size scaling ansaetze Eqs. (1) and (2) are of course never exact, and are usually understood as 

lim N-vp p=Pc . N (m = z/N*) = f{z) (1) 

N— >oo 

and 

lim (p-p c = z/N e )-l 3 (m) =g(z) (2) 

N— >oc 

for any fixed finite value of z. The limits here are pointwise limits, i.e. the norms of the differences between left and 
right hands converge uniformly to zero in any finite interval of z. Furthermore, f(z) and g(z) are usually analytic 
(holomorphic) for all finite z. 

In a typical first order (discontinuous) transition, an attempt to construct f(z) would give a function with two 
(5-peaks. In that case the convergence could at best be weak, i.e. in distribution sense. Usually one prefers to call this 
not finite size scaling at all, although this is strictly spoken a matter of taste and convention. 

In explosive percolation one has a "mixed" situation: For the right hand peaks in Figs. 1 and 2 of the main paper 
one has pointwise convergence, if one chooses rj = n + . But then the left hand peak converges to a <5-peak, i.e. the 
entire function f(z) is approached only in the weak sense. Similarly, for g(z) the convergence is strong (and g(z) is 
analytic) for z > (strict inequality!), where only the right hand peak of f(z) contributes. The function g(z) must 
vanish identically when only the left (<$-) peak contributes, which means that it must have a singularity at zq < 0, and 
convergence can only be weak in any interval containing zq. As for first order transitions, it is a matter of convention 
whether one calls this finite size scaling at all (as we did in this paper). 



ANALOGA TO FIG. 3 (MAIN PAPER) FOR THE OTHER THREE MODELS 

In the main paper, we showed in Fig. 3 for the da Costa model how (m) scales for p > p c , and we said that a data 
collapse similar to that shown in the inset holds also for the other two off-lattice models, while the collapse is much 
worse for the 2d model. We now show these data in Figs. SI to S3. 
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VARIOUS OTHER PLOTS FOR THE DA COSTA MODEL 



Although the scaling behavior of the order parameter with N at p — p c can, in principle, be inferred from Fig. 4 
(main paper), we show the data also explicitly in Fig. S4. In this figure we use three possible values of p c (one of 
them being the value obtained in 0], in order to show how strongly the exponent r]o depends on p c . 
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FIG. 3. (Color online) Same as Fig. SI, but for the AE model. This time the data collapse is better than in the other models 
(see the inset), reflecting the fact that the AE model is closest to an ordinary second order transition, among the four models 
studied here. 
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FIG. 4. (Color online) Log-log plot of (m) at p — p c plotted versus N. In order to show the sensitivity of the exponent r]o to 
the precise value of p c , curves for three values of the latter are shown. 



Results for p c (N), the effective critical points on finite systems, are given in Fig. S5. Notice that values of p c {N) 
depend crucially on the operational procedure used to define effective critical points. One possibility would be e.g. 
the point where the two peaks in P(s max ) have equal height (Fig. I). For the da Costa model, this would give 
non-monotonic dependence on TV. More natural seems the definition via equal areas under the two peaks. Notice that 
this would give ambiguous results for the 2d and AE models, as there the dips between the peaks are not very deep. 
But for the da Costa model this is unproblematic. Figure 5 shows that our best estimate of p c is slightly below the 
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value of [2|, giving thereby the largest contribution to the uncertainty of the exponent S (this slight inconsistency is 
also the reason why we used also these smaller p c values in Fig. S4) . 
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N 

FIG. 5. (Color online) Log-log plot of p c (N), defined as the value where the areas under both peaks in P(s max ) have equal 
height, versus N. Notice that there are only two points for N = 2 26 , since our p c (N = 2 26 ) is smaller than the p c value of Q]. 

Alternatively, we could define p c {N) as the point where P(s max ) has maximal variance. Variances of -P(s max ) /N 1 ^ 7 ^ 
are plotted in Fig. S6 against (p — PcjN 1 ^ 2 . We see distributions which are for small N markedly skewed and shifted 
away from the origin, but which become increasingly symmetric and centered at the origin as TV increases. Although 
this gives a much less precise estimate of S than Fig. S5, it demonstrates also that S > 0. 

Plots similar to Figs. S5 and S6 were not made for the other models, the main reason being the larger uncertainties 
of p c . 



FOR THE 2D MODEL, 6 IS STRICTLY SMALLER THAN 1/2 

As we said in the paper, one observation that corroborates = 1/2 for the off-lattice models is that it gives very 
"regular" behavior of (m) in the near subcritical region. Instead of showing here the evidence for this, we show for 
the 2d model what can go wrong, if is chosen badly. More precisely, we show in the top panel of Fig. S7 results for 
= 0.47, and results for = 0.5 in the lower panel. It seems clear that panel (b) is not very plausible, in particular 
since the scaling m_ ~ N~' n - with r)^- > rj + of the left hand peaks in Fig. 1 requires that the curves decrease with N 
for z < 0. For the other models similar curve crossings appeared when = 2 was used instead of = 1/2. 



FURTHER COMPARISONS WITH PREVIOUS PAPERS 

a) To our knowledge, the only previous work where 0^5 was seen in a continuous phase transition is Q. This 
dealt with interacting self avoiding walks in 4 dimensions. The authors also found double-peaked distributions of the 
order parameter, but they did not report different scaling laws for both peaks. Thus the reason for ^ S was not 
explained, as in our case, via a singularity of the scaling function g(z). Also, in 0] the authors found ^> 5, while 
we found in our models < S. 

b) A different scaling theory for the PR model was presented in [H ■ The authors there started from the assumption 
that (m) is independent of N at p = p c , which is definitely not true for any of the four models according to our 
simulations. Although this prevents our theories from being equivalent, there are some similarities. In particular, 
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FIG. 6. (Color online) Variances of s max for fixed N, plotted versus z — (p — p c )Af 1//2 , using the p c value of Q|. The width 
rescaling and the normalization are such that the curves should collapse for N — >• oo. For the N values shown in the figure, the 
collapse is far from perfect, but this is to be expected. Notice that the horizontal peak positions hardly change for N > 2 20 , 
showing that 8 > Q. 



the authors of Q show that the width of the FSS region scales as N~ e with 9 = 0.48, which is very close to our 
conjectured value = 1/2. 

c) A detailed study of the 2d model was made in The most remarkable agreement is that rjo was measured 
there as 0.0589(10), while we found 0.0612(8). The slight discrepancy is partially due to a slightly different estimate 
of p c (0.526565(5) in [H| against 0.526562(3) in the present work). Also the estimate of the Fisher exponent r in Q 
is fully compatible with our value of /3, if we accept the relationship between the two exponents given in 0]. 
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FIG. 7. (Color online) Linear-log plots of (s max )/A rl_/3e against (p — p c )N e for the 2d model. In both plots we used the values 
for P and p c given in Table 1. In the upper panel we used = 0.46 (the average between 0i and ©2), while we used = 0.5 
in the lower panel. 



